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1 INTRODUCTION 



Axisymmetric incompressible modes of the magneto-rotational instability (MRl) 
with a vertical wavenumber are exact solutions of the non-linear local equations of 
motion for a disk (shearing box) . They are referred to as "channel solutions" . Here, 
we generalize a class of these solutions to include energy losses, viscous, and resistive 
effects. In the limit of zero shear, we recover the result that torsional Alfven waves are 
exact solutions of the non-linear equations. Our method allows the extension of these 
solutions into the dissipative regime. 

These new solutions serve as benchmarks for simulations including dissipation 
and energy loss, and to calibrate numerical viscosity and resistivity in the Zeus3D 
code. We quantify the anisotropy of numerical dissipation and compute its scaling 
with time and space resolution. We find a strong dependence of the dissipation on 
the mean magnetic field that may affect the saturation state of the MRl as computed 
with ZeusSD. It is also shown that elongated grid cells generally preclude isotropic 
dissipation and that a Courant time step smaller than that which is commonly used 
should be taken to avoid spurious anti-diffusion of magnetic field. 
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The magneto-rotational instability (MRl) l|Balbus fc Hawlevlll991h . is generally regarded as the best candidate for explaining 
the "anomalous visc osity" that transports a ngular momentum through accretion disks. The first two-dimensional numerical 
studies of the MRl (jHawlev fc Balbuslll992l ) revealed a surprise: the fiow broke along the vertical axis into two distinct and 
regular sheets, which appeared as channels when visualized in a meridional plane. In this plane, which projects out azimuthal 
motion, the perturbed velocities are nearly radial, one fiowing inward, the other outward. This pattern was given the name 
of channel solution. It is a recurring feature of more general three-dimensional global numerical simulations of the MRl, 
appearing intermittently. Indeed the most unstable linear modes of the MRl are modes with a vertical wavenumber, which 
are, in effect, channel solutions. 



Goodman fc Xu 1 19941 ) showed that these modes could be destroyed by three-dimensional parasitic instabilities, notably 



the magnetic Kelvin-Helmoltz (KH) instability. Such studies are important for the light they may shed on the process of the 
saturation of MRl turbulence. Resistivity and viscosity are also key to the saturation of the MRl, as channel flows break 
down and fields reconnect. As the power of computers increases, more of these important microphysical processes are being 
included in the simulations. On the other hand, as the complexity of codes increases, fewer benchmarks are available to check 
and calibrate the implementa tion of numerical metho ds. It is this issue that motivates this paper. 

One of the key insights of lGoodman fc Xul (|l994r ) was to note that the finear modes of the MRl are exact solutions of the 
non-finear local equations of incompressible magnetohydrodynamics (MHD). This is also seen in non-linear circularly polarised 
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torsional Alfven waves, whose linear modes exactly satisfy the non-linear MHD equations for a polytropic (or incompressible) 
gas. In this paper, we will establish a connection between these two types of solution. 

It is the purpose of this paper to extend the Goodman-Xu solutions to include viscosity and resistivity. We are also able 
to find solutions that include a net heating or cooling term, provided that such gains/losses are linear in the temperature. We 
first compute the linear MRI modes with a vertical wavenumber in the shearing box regime with viscous and resistive terms 
included. We show that the only remaining non-linear term is the total pressure gradient (magnetic plus kinetic), and we 
establish a condition for it to vanish (Section 2). In the next two sections we investigate less general cases without viscosity 
(Section 3) or without shear (Section 4), for which we can give extended families of analytic solutions. In Section 5 we provide 
methods for finding isolated solutions under more general assumptions. We then present two applications of our results: in 
Section 6 we benchmark a new version of the ZeusSD code with a conservative scheme for total energy; in Section 7 we use 
our solutions to compute the numerical resistivity and the numerical viscosity in the ZeusSD code. We discuss our results and 
conclude in Sections 8 and 9. 



2 GENERAL METHOD 
2.1 Shearing sheet 

The shearing sheet system results from a local first order expansion of the dynamical equations of m otion, with the inverse 
radius serving as the small parameter. This approach dates back to a celestial mechanics calculation of Hill 1 18781 ). The frame 
of reference rotates at circular angular velocity fl. Radial, azimuthal and vertical directions are labeled by local Cartesian 
coordinates x, y and z. The origin of the frame follows an unperturbed fluid element moving in a circular orbit. The radial 
logarithmic derivative of is 

2dhrR'"=° 

and characterises the local shear. 

The fundamental dynamical equations in this rotating frame are the mass continuity equation, 

g + Vipv) = 0, (1) 
where p is the mass density of the gas and v is its velocity, and the Navier-Stokes equation with a kinematic viscosity vv, 

^ + {vV)v + 2nzxv + V{2AQ.x^) + i V(p + ^) - -(B-V)B = -V-(pvvcr) (2) 
ot p 2 p p 

where p is the thermal pressure, B is the magnetic field (divided by 2-^7?) and aij = ^{diVj + djVi) — \dkVk5ij is the stress 
tensor. The vertical gravity is neglected. The induction equation is 

^ = Vx{vxB -t^bJ) (3) 

where J = ^ X B and r^B is the resistivity. Finally, we adopt an ideal gas equation of state with adiabatic index 7 so that the 
internal energy equation reads 

^ Dln^ = ,bJ' + puya: Vv - A (4) 

where A is the net cooling function. 

Equations ([T]) to (|4]) form the governing system of which we seek particular solutions. 



2.2 Incompressible vertical modes 

We seek solutions that are a single Fourier mode with vertical wavenumber k and a constant vertical magnetic field Bo 
superimposed on the background shear. The solutions have the form 

V = 2Axy + u B = Bqz + b (5) 
with 

r st-\-ikz 1 ri stA-ikz fn\ 

u — due b = dbe (6) 

We assume 5w _L i and 5b _L z. It is understood that the physical solutions are obtained by taking the real part of 
equations ((5)1. (In particular, note that the S amplitudes may be complex. ) 
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For this form of solution, the mass continuity equation states that the Lagrangian derivative of p is zero, hence the 
density remains constant along the trajectories of the fluid elements. Thus, p is constant in time and space provided that it 
was uniform initially, which we shall assume. 

We further assume that the pressure p is initially a function of z only. This property is also conserved in time for our 
particular flow. Upon substitution of equations ([5} and ([6]) into equation ((2|, we obtain 

su + 2Auxy + X u+-zdz(p + _ -Boikb = -uyk^u (7) 

P \ ^ } P 

where we assume that vy is uniform and we use the symbols SRfz] and to denote the real and imaginary parts of a complex 
number z. Note that this equation makes sense only if the z pressure gradient terms vanish, a key point to which we shall 
return below. 

Finally, assuming a uniform resistivity 773, the induction equation becomes 
sh — 2Abxy ~- Boiku = —rjBk^b (8) 
In what follows, we shall use the new variables 

; 2 

u — k uv 
and 

; 2 

1] — k TIB- 

The only remaining non-linear term in the above equations is the total pressure gradient dzPtot with ptot ~ p+ (5R[b])^/2. 
It is also the only term of this equation with a z component. If, for the moment, we assume that this term is zero, then the 
problem for u and b is contained within equations ([7]) and ([8]) and decouples from the energy equation. We are left with a 
simple linear, incompressible problem with viscosity and resistivity, which we now turn to solve. 

The linear Euler equation may be written 

Eu^^b (9) 
P 

with the 2x2 matrix E defined by 

E = ^,1 ' (10) 



s + u -20, 
2(n + A) s + iy 



where we considered only the x and y components of the vectors. The induction equation can be rewritten similarly 

¥b = Boiku (11) 

with 



F = 



s + i] 
~2A s + r/ 



(12) 



Operating on equation ((TTJ with E and using equation ® produces the linear eigenvalue problem 

{E¥ + k'^VAl)b = 0, (13) 

where I is the identity matrix and k'^v\ = B^k^ / p where va is the Alfven speed, s is therefore a root of the determinant (det) 
polynomial 

P^,,,(s) = det[EF + k'^v\l], 

which is precisely the MRI dispersion relation. 

We need compute directly only the restricted case Po,i7(s) for u = since the full linear problem with parameters (s, v, rf) 
is equivalent to the one with [s + v,Q,rj ~ v). Hence, P^,n{s) — Po,i7-i/(s + i^)'- 

Po,^(s) = (r; + s)'(k' + s^) + 2 [2An + s{s + r])] k^v\ + {kvA^ = (14) 
and in developed form 

Po,„(s) = s* + 2773^ + (77^ + 2k'^v\ + k)s^ + 2r){k'^v\ + k)s + kt)^ + 4Anfe%i + {kvAf (15) 

with — AQ.{A + Jl) which is identical to the form given by equation (12) in [Fleming et Zl (|2000l ). The general dispersion 
relation may then be obtained by replacing s hj s + u and r; by 7; — in equation (I14|l : 

P.,r,{s) = (77 + sf [k' + {u + sf] + 2 [2An + (s + u)(s + 77)] fc^4 + [kvA^ = (16) 
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2.3 Condition for a homogeneous total pressure 

For the solutions whose form is of the previous section, the internal energy equation becomes 

^ = (7 - 1) [vmb]f + pumu]f - A] . (17) 

Note that the v-'Vp term vanishes because has components only along x and y whereas Vp is along z. The above may 
be rewritten as an equation for the total gas plus magnetic pressure: 

^ (p + = ?fi[sb]-m + (7 - l)M^[b]f + piy{Q[u]f - A). (18) 

From here on, we restrict the cooling function to be of the form A = —T + ap with F and a constant. Since p is a constant, 
this is equivalent to taking A to be a linear function of temperature. (In effect, this is just the leading Taylor series expansion 
of A(r) around an arbitrary point in temperature.) The effect of F may be absorbed into ptot by adding a linear function of 
time with no spatial dependence, and without loss of generality, we may set F = 0. We rewrite the time evolution equation 
for the total pressure accordingly: 

^ + (7 - l)aptot = 5R[s6]-K[b] + (7 - l){rimb]f + piymu]f + |K[6].K[b]). (19) 

The right hand side of equation (|19|l must be independent of position for a self-consistent solution. Formally, these terms may 
be expressed as a spatially constant term plus a term of the form ^[aexp{2st + 2ikz)], with the complex amplitude a constant 
in both time and space: 

a= [s-(7-l)(77--)J ^^{j-l)up ^ (20) 

where we remind the reader that Sb^, Sby, Su^ and Suy are complex numbers. Our task, therefore, is to investigate the 
conditions under which a vanishes. If this can be done, the full set of equations is reduced to the linear problem of the 
previous section plus equation p9p . which now reads 

+ (7 - l)«Ptot = + (7 - l){rj + -)J 2 + ^'^ ~ ^'^^■^ 2 ^ ^ 

where \Z\^ — (5R[Z])^ + (9[Z])^ (modulus). When the solution of the linear problem ((Tj-© (without the total pressure 
gradient) is inserted into this last equation, one gets a simple solution for ptot of the form : 

Ptot(t) = Pi exp(2K[s]t) + p2 exp(-(7 - l)at) (22) 
with 

" (7 - l)a + 2K[s] 
and 



[m + (7 - + f )] ^M^^ + (, _ D.p i'^--! (23) 



P2 =Ptot(0) -PI. (24) 

The above is a solution of the non-linear problem if and only if the total pressure gradient term vanishes, hence if and 
only if a = 0. In the case with shear, this condition can be recast in the form of a fifth order polynomial Q{s) = 0. To see 
this, we first express the components of 5b and Su in terms of Sb^ only. The first row of equation (I13|) gives Sby as a function 
of Sb^c'- 

,Sby_^ {s + u){s + r])+iAn + k\\ 
~ 5b^ 2Q{s + ri) ' 

Note this is valid only provided that Q is non-zero. If fl — A = 0, then Sbx and Sby are independent (see the case without 
shear in section |4] below) . 

Equation (|ll|l now allows us to express Su in terms of Sb and hence Sb^: 

Sux = -rr^{s + ri)Sbx (26) 
and 

= ttV (s + ry)/]5b,. (27) 

a can now be rewritten 



a = Sbi 



[s - (7 - l)(r? - |)] i(l + f) + Sbl (7 - 1)^^^ [{s + vf + (-2A + (.s- + r,)ff] . (28) 
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Finally we substitute / thanks to equation psp and gather quantities over the same denominator: 



7-1 Q{s) 



2 (s + r?)2 
where Q{s) is the following polynomial in s 

1 



(29) 



Q(s) = 



1 a 



+ 



(s + V) + 



{{s + u){s + 1]) + k\l)^ 



(30) 



The homogeneity condition a = is therefore satisfied if s is a root of Q. Henceforth, we refer to Q as the homogeneity 
polynomial. Note that it is generally of order 5. 

The parameters of the system are hence Q and A for the shear, kvA (which actually combines So, k and p) for the 
magnetic field, i/, rj, 7 and a (or A) for the properties of the gas. For a given set of these parameters, we now want to find a 
growth rate s which satisfies the dispersion relation (|16|) and for which the total pressure gradient vanishes, ie: P{s) — and 
Q(s) = 0. 

We shall first restrict our analysis to real roots s in simple cases, although we treat the case with no shear exhaustively 
(see section |4ll . Most of the solutions presented are therefore standing wave solutions, except in the case without shear. There 
we find circularly polarised waves and other propagating disturbances. 



3 INVISCID SOLUTIONS 

In this section we set v — 0. The homogeneity polynomial becomes 

(31) 



We choose to find a common real root to P and Q. The only way Q can have a real root is if s = (7 — 1)(?7 + ^) because the 
quantity in the square brackets of expression (|3ip is strictly positive for s real. 

3.1 Cooling/heating 

The extra degree of freedom granted by the presence of thermal losses/gains makes such solutions easier to find: given a real 
root s of P (ie: a standing mode), we need only adjust a to 

a = 2(,7--^s). (32) 
The final solution is then simply given by the expressions (jSj, (O and (|22p . 

3.2 a = 



From (|32l) with a = we immediately see that s = (7 — 1)77 is required for a uniform total pressure. Using this in the relation 
P{s) — 0, we obtain a quadratic equation for 77^: 

7''(7 - 1)%^ + 27[(7 - i)k^yl + 2n7(^^ + A)]rj'' + k'^vliAAfl + fc^^i) = (33) 
Alternatively, this can also be viewed as a quadratic equation for k'^v'\: 

{kvAf + [^AVl + 2(7 - l)-irf]k'^v\ + -I'^rf [m{A + f7) + (7 - l)2r?2] = 0. (34) 

Either of these equations allows a determination of the set of parameters {ri,kvA) for which there exists a solution. For 
example, in the case where = + ^) > 0, equation (|34|l has a real root for k'^v\ provided that 

2 QA^ 

which sets an upper limit on the resistivity. The request that this root be positive sets up the additional constraint 

2 -2An 

which forces A > and sets an additional upper limit on the resistivity. The final condition (an upper limit on resistivity) is 
2 _ -An . / 2 -A 



V < min -, . (37) 

7 V7-l^ + 7^^. 
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3.3 rj = V — 0, adiabatic 

In this case, the condition p4|l simply becomes kvA = 2\/— Afi, so ^4 < for such a solution to exist. The growth rate 
is then s = 0, a marginally stable mode of the MRI. The dispersion relation has only real roots, s — (double root) and 
s — ±2-\/r2(f2 — A), but the growing and decaying modes do not fulfill the homogeneity condition. (As mentioned above, this 
solution is also valid when the constant F is non zero, which is not, strictly speaking, adiabatic.) 



4 NON ROTATING FLOW 



Here we set A = Q, = and drop the assumption that s is real. Without rotation and shear, x is no longer a special direction; 
the direction z is still defined by the mean field Bq. The system is now invariant under rotation of axis z and the eigenvectors 
of the linear system now depend on two independent variables, say Sbx and Sby. The effective dispersion relation becomes: 

\2 



P(s) = R(sy = 
with 

R{s) = (s + 7]){s + + k^v\. 

Without shear, the homogeneity condition a = with the definition (|20p becomes 



7- 



1 



7-1 



+ 



Sbl + Sbl 



- 



(38) 



(39) 



(40) 



This can be achieved if either Sb^ + Sby — 0, or if the factor inside the brackets vanishes. 



4.1 Torsional Alfven waves 

In this case, we simply set 5b^ + Sby = by choosing either of the circularly polarised cases Sb^ — -iziSby. Now, both roots s 
of the dispersion relation R = Q provide a possible solution: 

^±--^±y(^)^-fcM (41) 

When kvA < \r] — v\l2 we get two standing modes. When kvA > |?? — i'|/2, s± have imaginary parts and the two 
solutions correspond to right or left circularly polarized waves. In particular, when a = — rj — 0, we recover circularly 
polarised torsional Alfven waves which are indeed well-known solutions of their non-linear governing equations. 



4.2 Non polarised waves without shear 

If Sb^ + Sby is non zero, then we need to find the common roots of P and the simple quadratic 

Q(.) = f-. + :^ + ^(3 + .)^ (42) 

which is the factor inside the brackets of equation (|40|l . Since P = R', to find a common root of P and Q means to find a 
common root of Q and R. For simplicity we assume 7 — 5/3, but it is not much more difficult to do without this assumption. 

We detail our analysis of the common roots of P and R in appendices |A] for complex roots and[B]for real roots. Here, we 
simply summarise our result that a > y/6\kvA\ is a necessary and sufficient condition for the existence of common complex 
roots, i.e. : the existence of propagating disturbances as solutions. We are also able to give an expression for kvA in terms of 
the other parameters of the problem in the case when there exist a common real root, i.e. : when a standing wave is solution 
of the problem. 



5 SOLUTIONS WITH SHEAR, RESISTIVITY, VISCOSITY AND COOLING 

In general, one may not be interested in the complete range of parameters for which a solution exists. A benchmark calculation 
only needs one set of parameters. In that case, we may simply pick a growth rate and treat P = and Q = as equations 
for k'^VA (both quadratic). Then the process of finding a set of parameters that yields a solution is greatly simplified. As an 
illustration, we set Q, = 1, A = —3/2, u = 1/5 and r] — 1/10 and seek k'^VA and a as functions of s. P = implies 

fc^«i = ^ (74 - 5s(3 -f 10s) ± 5^218 - 10s(ll + 40s)) . (43) 
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Then the equation Q = is linear in the variable a. For example, if we now seek a standing wave solution with the growth 
rate s = 1/2, we find 

= ^i±l^ . 1.874 (44) 
and from Q(l/2) = 

899569 + 13560^/7 ^ 
597490 

Here we have an example of an explicit benchmark with viscosity, resistivity and heating in a shearing box. 



6 NUMERICAL BENCHMARKS 

The original ZeusSD code (see lstone fc Norman 1992al lbh is not written in a fully conservative form. In particular, equation 



Q is used to compute the evolution of internal energy. In general, this scheme leads to significant loss of total energy. For 
example, if discretisation errors lead to kinetic or magnetic energy losses, this artificial dissipation is not reflected in viscous 
heating or Ohmic heating and the total energy decreases. This energy is effectively "radiated" away. 

However, when viscous and/or resistive terms are included in the code, part or all of the total energy loss is recovered as 
heat and the total energy loss is reduced. As an illustration, we ran torsional Alfven tests with the original ZeusSD code, and 
with a total energy conserving scheme (see below). Fig.[T^ shows that the total energy scheme performs much better in a case 
wi thout explicit resistive and viscous terms: the internal energy scheme heavily distorts the wave profile (this was already noted 



by Turner et al.| [2003). On the other hand, fig. [T]d shows that the internal energy scheme with some resistivity and viscosity 



is indistinguishable from the total energy scheme. Because of the finite resolution (32 zones) there is still some numerical 
dissipation and the numerical results are slightly damped compared to the analytical solution. The difference between the 
analytical solution and the actual simulations disappears on the scale of these graphs at a resolution of 128 zones for both 
these tests. 

In principle, grid based schemes cannot avoid numerical dissipation. However, it is possible to make numerical dissipation 
look more like physical dissipation by using a conservative form for the evolution equations. For example, we can evolve the 
total energy, and deduce the internal energy as the remainder of the mechanical (kinetic, magnetic plus potential) energy 
subtracted from the total energy. In this case, we write the total energy equation as: 

f + V.^=-A (46) 
with 

£ = e + ^pv^ + + p$ (47) 
where $ — 2AQx^ is the tidal potential energy, e is the internal energy and the total energy fiux is 

1 + e + ipi;^ + p<i>j + {Bxv)xB + tibJxB - pwcr-v. (48) 



We have implemented this in the Zeus3D code. This is similar to the work of lTurner et al. I l|2003l) andlHirose et al.l (120061) 



but we also include Ohmic stresses rjsJxB and tidal potential energy fiux vp$. The question arises at what stage of the 
calculation each term of the total energy fiux should be evaluated. We ran various benchmarks (torsional Alfven waves 
and standing mode solutions presented in the previous sections) and varied the order with which the fiuxes were computed. 
We noted that it is crucial to compute each fiux term simultaneously with its corresponding source or transport term. In 
particular, it is critical to compute (Bxv)xB us ing the time centred values for B and v computed with the method of 



characteristics (MOC, see lStone fc NormanI 1 1992131 ). On the other hand, the kinetic energy fiux should not be directionally 
split the way the momentum transport step is. 

Finally, it is much better to join the tidal fiux to the density transport term and not to the tidal force source term. We 
now illustrate how we used our analytical solutions to prove this last point (see fig. [2| . We used our code with two slightly 
different versions in order to reproduce the analytical solution presented in section [5] The first version (dotted lines on fig.[2ll 
would compute the tidal potential fiux p^v at the same time as the tidal source term. The second version (dashed lines on 
fig. [2]) would compute this fiux jointly with the transport step. In the first version, the resulting total pressure gradient is not 
fiat and the magnetic energy loses its low z/high z symmetry. The second version retains the correct symmetry and displays 
a fiat pressure profile. Note however that in both computations the average total pressure and the magnetic pressure are 
slightly lower than the analytical solution. Both simulations shown are for cubic boxes of 32 zones aside and higher resolution 
improves the magnetic pressure more efficiently than the total pressure. 

The final scheme we adopted was to compute in five distinct steps: 



8 Pierre Lesaffre & Steven A. Balbus 

(a) Alfven Test r]g=Uy=0 (b) Alfven Test r]g=0.Q5 Vy=0.lQ 




Z Coordinate Z Coordinate 

Figure 1. Torsional Aflven waves tests. Solid line: analytical solution, dotted line: internal energy scheme, dashed line: total energy 
scheme. Parameters are fl = A = a = 0, -f = 5/3, So = = 10, k = 2it. The simulation box is a cube of 32 zones aside and has 
physical length 1. The plots show snapshots of the azimuthal magnetic field evaluated on a vertical line, (a) Left panel: i/y = r]g = 0, 
after 3 oscillations, (b) Right panel: vy =0.1 and -qs = 0.05 at a time correponding to 9.9 oscillations. 



(a) Tidal flux position test 

6201: 



610 r 

(u - 
^ 600 - 

cn 
CO 

CU : 



o 580 T 

H 

570 I- 



5R0 P ... I ... I ... I ... I 
0.0 0.2 0.4 0.6 0.8 

Z Coordinate 



(b) Tidal flux position test 
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Figure 2. A benchmark with shear, resistivity, viscosity and heating (see section[5]l. The wave number is fc = 27r and p = 1. We display 
the total pressure (a, left panel) and the magnetic pressure (fe^ + by)/2 (b, right panel). The time is i = 6 which corresponds to 12 
e-folding times (s = 1/2). The solid line is the analytical solution, the dashed line is for a total energy scheme where the flux of tidal 
potential p^v is computed jointly with the transport step and the dotted line is a total energy scheme where this flux is computed with 
the tidal source term. 



• vp is first computed using an upwinded pressure computed at the same time as the pressure gradient source, 

• the viscous term is computed at the same time as the viscous forces, 

• the remainder of the flux w (e + \pv'^ + p<E>) is added after the hydrodynamical transport term, 

• the resistive term is computed along with the resistive electromotive force, 

• the {Bxv)xB term is finally computed with the MOC advanced v and B which are used for the constrained transport 
of B. 

Along this process, we evolve the internal energy e thanks to equation Q. In particular, this provides an advanced 
estimate for e in the flux term v{e +...). At the end of these steps, we compute £* with the updated values of all variables. 
We then use equation (|46p with A = to advance the total energy to its new value £. If the internal energy scheme was 
perfect, we would have £* = £. However, this is almost never the case and a correction £ — £* needs to be applied to e in 
order to conserve total energy. We deduce the rate of correction of internal energy e — {£ ~ £*)/-U where At is the length of 
the time step. The internal energy is flnally updated with A 7^ and e thanks to an isochore heating/cooling step. 
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7 NUMERICAL VISCOSITY AND RESISTIVITY 

In appendix [C] we present a method to estimate the numerical resistivity and viscosity in a code. The idea is to probe the 
numerical dissipation in the absence of explicit resistivity and viscosity, and to determine the effective numerical dissipation 
coefficients by fitting results to our analytical solutions that include viscosity and resistivity. In the following, we write for 
short ?;n and for the numerical resistivity and viscosity. As explained in the appendix[C] we measure directly (r;N + VN)k^ 
and {jiN — UN)k'^ respectively to second and first order in (t^n — VN)k/vA- We then deduce the values for j^n and i'n. We used 
this method mainly on the internal energy scheme version of the ZeusSD code since it is the version that is generally used in 
published applications. 



7.1 General trends 

Our method provides a direct estimate for the numerical dissipation in a code. Therefore it gives the numerical floor for the 
physical viscosity and resistivity in a given code. For codes devoid of a viscous or resistive term, it also allows to compute the 
effective Reynolds and Prandtl numbers. We now investigate general trends of the numerical dissipation. 

We first examine wave numbers along the vertical direction. In figure [S] we examine the dependence of t^n and !/n with 
various parameters. The Courant number (or Courant coefficient) is a parameter that controls the time step of a code. In the 
ZeusSD code it is defined as 



max v'^ + + vYj (49) 

where v, c and va are the local speed, sound speed and Alfven speed in the fluid, Ax is the size of a zone. At is the size 
of a time step and the maximum is taken over all grid zones. We measured the dependence on resolution, wave number, 
perturbation amplitude and mean field amplitude for three different Courant numbers: 0.01, 0.1 and 0.5. we display the 
results only for a Courant number of 0.1 and we discuss the differences when applicable. We first ran a standard run with 
parameters f3 — 2/Bq — 400, k = 2-k, an amplitude of \Sbx\ = 0.001, p = 1, a Courant coefficient of C = 0.1 and a spatial 
resolution of 32 zones in all three directions (hence Ax = 1/32 since we use a physical length of 1 for the size of the box). We 
then varied each parameter in turn away from these values. 

Figure [3^ and [Sja show that rj-^ and scale linearly with the size of the time step and as the square of the size of a 
grid cell. An interpretation of these trends is that our scheme is 1st order in time but 2nd order in space. Note that at a 
Courant coefficient of 0.5, the numerical t^n changes sign. As a whole, the numerical scheme remains stable in the sense that 
?7n + i^N is always positive. However, rjt^ or individually could be negative. ?7n < indi cates that the MHD part of the 
time step behaves like antidiffusion. Antidiffusivity in Zeus_was already noted bvlFallj l|2002l ) who also pointed out that lower 



Courant numbers lower antidiffusion. More recently, iFromang &: Papaloizoul (|2007l ) also pointed out antidiffusion in Zeus at 



large scales. Here, we quantify the effect in more detail. The wave number with the lowest numerical resistivity turns out to 
be fc = 2-k{x + y + z). The resistivity of this mode is negative for all Courant coefficients above 0.12 (see dashed line on figure 
[Sjj). Such negative values for the resistivity are only found for wave numbers with coordinates lower or equal than 2: only the 
largest scales are affected. With the Zeus3D code, it might nevertheless be safer to adopt Courant coefficients below 0.5 or to 
include some minimal amount of physical resistivity in the code. Including physical dissipation has the advantage that it will 
also improve the energy budget, as noted in the previous section. 

In figure ^jp) it appears that the dissipation has a finite limit as the time step tends toward zero. Indeed, the finite space 
resolution does not allow the scheme to achieve an infinite precision. Similarly, in figure ([3^) the scaling of the numerical 
dissipation is a power of -2 in the number of zones at low space resolution, but at high space resolution it turns into a shallower 
power of -1. Indeed, the Courant number is kept fixed (hence At/ Ax is fixed) and the scheme is second order in space but 
only first order in time: at high resolution, the numerical dissipation is dominated by the order of the time integration scheme. 
At a higher Courant number of 0.5, the shallower slope of -1 occurs at even lower space resolution. At a Courant number of 
0.01, the slope of -2 is seen over the whole range of space resolution we tested. Note that the numerical resistivity also turns 
out to be negative for the highest resolutions at Courant numbers 0.1 and 0.5. 

Figure [Sj; shows the dependence of t^n and v-^ on the wave number. Unlike a physical viscosity or resistivity, t^n and i^n 
vary according to the length scale, with a maximum at fe = 8 x 27r£ (four grid points inside each wavelength). There is no 
clear scaling but figure |3}; suggests that the total dissipation behaves roughly like a power with an exponent between 1 and 
2 (1.6 seems to be the best match) until it reaches the maximum dissipation. The quantity (t^n — i'N)k/vA can be as high as 
0.5 for k = 8 X 2-kz, so the method of appendix[C]is not accurate for higher wave numbers. 

Our estimates show that immcrical dissipation is nearly independent of the amplitude of the initial perturbation (although 
there is a slight dependence on it at a Courant number of 0.5). However, we observe a strong dependence on the amplitude of 
the mean magnetic field (see figure[3ji). The numerical viscosity appears to be directly proportional to the mean magnetic field 
at low Courant number (0.01) with some additional dissipation at low {3 for higher Courant numbers. The numerical resistivity 
obeys the same law, with an additional change of sign at low values of 13 (note: it remains positive at low (3 for a Courant 
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(a) Dependence on Space Resolution 




(b) Dependence on Time Resolution 
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Figure 3. We plot the numerical viscosity (solid) and resistivity 77N (dotted) vs various param eters such as (a) total zone number 
(the dashed line is tjn + for fc = 27r X 8z), (b) Courant coefficient fsee lStone &: Norman|[l992al lbl)(the dashed line is r)^ for the wave 
vector fe = 2tt{x + y + z) and the dash-dotted line indicates the zero threshold.), (c) kz for a mode with k along z and (d) amplitude of 
the mean field (the dashed line is tjn + i^N for fc = 27r X 8z) 



number of 0.01). This explains why the benchmarks of the previous section (which are for large mean field amplitudes) have 
strong dissipation compared to the standard runs of the present section. This might be an issue for the computation of the 
saturated state of the MRI with ZeusSD. If numerical dissipation does increase with the turbulent magnetic energy, this could 
affect the total effective dissipation in the system. 

To summarise these results, we suggest to approximate the total numerical dissipation with the following scaling formula; 

/ , s 1.6 

7?N + i^N^0.76Aa;^/3"M _j + 1.08Aa;C/3"\ (50) 

We calibrated both coefficients of this formula on figure |3Jd and the exponents for Aa;, C, k and (3 are obtained from figures 
EKi Et, |3j; and[3}l respectively. Formula (|50p should therefore be taken only as indicative for values of parameters not too far 
from those tested here. Furthermore, as stated, the scaling in k should also be taken with caution (see figure |3j;). 



7.2 Anisotropy 

Our method allows us to quantify the anisotropy of the numerical dissipation. We measured the numerical dissipation for all 
wave vectors in the Fourier domain of the box with coordinates of the form kt = 27rn with i = x,y,z and ^ n 16 (a 
grid of 17'' — 1 measurements). Many shearing box simulations actually use half the resolution in the azimuthal y direction 
compared to the radial x and vertical z directions. We therefore did the same measurements (with fcy ^ 8 x 2-k) on a cubic 
box with 32x16x32 zones in which the grid cells have an aspect ratio 1:2:1. 
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7.2.1 Cubic grid cells 

In figure [4^ we plot all t^n + i^n measurements for the 32'^ (cubic cells) simulation against the norm of the wave vector. The 
overall shape of this diagram roughly follows figure [3]; with a maximum of dissipation at 15x27r. Even for cubic grid cells, 
the numerical dissipation already shows some degree of anisotropy: at a given wave number it varies widely. We detail the 
distribution of this spread in figure [5K for wavenumbers k which have 2nNk ^ k < 2n{Nk + 1) with A^fe = 15. Wave vectors 
with the highest dissipation are those that point towards a cartesian axis. For a fixed jfcj, wave vectors along an axis maximise 
the size of a single component. We therefore suggest that the numerical dissipation at a given wave vector is dominated by 
the dissipation at its maximum coordinate. For wave vectors of norm \k\ higher than 16, the three coordinates have similar 
values, hence the numerical dissipation is more and more isotropic. Interestingly, the numerical Prandtl number 

PniN — m/riN 

is quite isotropic for all wave numbers and slightly decreases from 2 at small wave numbers to 1 at large wave numbers (see 
figure [l]:; the very small wave numbers have higher Prandtl numbers, but the numerical dissipation is much lower there). The 
isotropy of the Prandtl number is even better at lower Courant coefficients (C = 0.01), with a mean value closer to (slightly 
above) 1 and a spread between 1 and 2. It is interesting to compare these results to the recent work of iFromang et al.l (|2007l ) 
who estimate Prandtl numbers between 2 and 4 for Zeus. It is also striking that all our measured Prandtl numbers are greater 
than 1. 

As mentioned, a few directions yield a negative resistivity. The corresponding wave vectors at C = 0.5 have their 
coordinates amongst the following list: (1,1,1), (1,1,0), (2,1,1), (2,2,1), (2,2,2) and their permutations. 



7.2.2 Elongated grid cells 

For elongated cells, the diagram |41d is only slightly more complicated. It is similar to figure UK, but replicates its pattern 
extended by a factor 2 in amplitude and squeezed by a factor 2 in wave numbers. This additional feature results from the 
halved resolution in the y direction. As shown on figure[5j3, the dissipation is not symmetric to x-y exchange. On the contrary, 
y wave vectors undergo much larger dissipation. This shows up even more at smaller wave numbers as seen in figure [5]; and 
[5ll. However, the Prandtl number does not show more anisotropy than in the case of a cubic cell: numerical resistivity and 
viscosity react in the same way to the resolution loss in the y direction. 



7.3 Scheme 



We tested our code in various configurations to investigate the impact on numerical dissipation. We found that isothermal 
simulations are slightly less dissipative than adiabatic simulations (with numerical resistivity more negative in general). We 
could hardly see any difference between the internal and the total energy schemes. We also found that to use the non-linear 
artificial resistivity as coded in IStone fc NormanI (|l992a[ l did not change the numerical dissipation: our tests did not trigger 
significant artificial viscosity because of the incompressible nature of our test fiows. 



8 DISCUSSION 

8 . 1 Incompressibility 

We discuss here a few caveats, limitations and possible extensions of our method. First, we wish to stress that the condition of 
incompressibility on the modes is a crucial one. Indeed, any change in density will alter the 1/p factor in the Euler equation, 
introducing additional non-linearities that might be difficult to address with analytical tools. 

In section 7 we have measured an equivalent numerical viscosity for torsional Alfven waves only. This gives a first estimate 
of the numerical dissipation, but an arbitrary MHD flow cannot be decomposed into such modes. For example we cannot 
probe any viscosity associated to compressible flows (see section . On the other hand, the viscosity of compressible fiows 
could be measured by studying the width of shock fronts, or with damped magnetosonic waves in the linear regime. 



8.2 Boundary conditions 

Periodic boundary conditions in z are essential for our analysis. For example, reflective boundary conditions will mix two 
Fourier modes which very likely will open the way for a cascade at many other wave numbers. 

It should be noted that our analytic shear solutions do not strongly test the implementation of shearing box boundary 
conditions. Indeed, except for the mean steady flow, our solutions for non-zero Q depend only on the z coordinate. For example, 
these benchmarks could not tell if the code uses periodic boundary conditions in the x direction or shearing box boundary 
conditions. 
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Figure 4. We show the total numerical dissipation tjn + '^N (panels a and b, upper side) and numerical Prandtl number Prnfj = uj^/rn^ 
(panels c and d, lower side) vs. k for every wavevector of the computational box. Panels a and c (left hand side) are for a 32x32x32 zones 
computational box. Panels b and d (right hand side) are for a 32x16x32 zones box with elongated cells. 



To improve this, we could need to find solutions with spatial variation in more than one direction. A non-zero kx is in fact 
perfectly tractable, and the equations hardly change if one uses the expression [k.Bo)'^ / p — k^B^/p instead of k^v\. However, 
it will not probe more efficiently the shearing box boundary conditions: periodicity in x would still be indistinguishable from 
shearing box conditions for most variables. 

In order to probe the shearing box boundary condition s, a non-zero azimuthal wave number is needed. Unfortunately, a 
non-zero ky yields an Eulerian wave vector changing in time IjBalbus fc Hawlevlll992r i. In that case, the homogeneity condition 
changes in time and the total pressure gradient cannot be dealt with at all times. However, it is worth noting that the total 
pressure term is still in this case the only non-linear term. Semi-analytic solutions of the equations without the pressure term 
can be found for MHD shearing waves. We plan in future work to benchmark the MHD shearing box boundary conditions by 
using such equations. 



8.3 Thermal diffusion 

A thermal diffusion coefficient can very easily be included in our analytical solutions. Under the assumption of a uniform 
density, the thermal diffusion term in equation (|19p is proportional to 

xAp = xA {ptot - ]ft[b]-^[b]^ = x'ik^m-m (51) 

where x the uniform thermal diffusion coefficient. We recall that the term due to cooling in equation (|19p is a/2 K[6]-5R[b]. 
To include thermal diffusion in our formalism is hence equivalent to use a -\- ik^x in place of a. 
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(a) rjf^+Vf^, 32x32x32 simulations, N^=15 
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(d) rjf^, 32x16x32 simulations, N^=5 




Figure 5. We show the total numerical dissipation rj-si + ussi (panels a and b, upper side) and numerical resistivity uyi (panels c and d, 
lower side) for each wave number with positive coordinates and such that 2ttNi^ fe < 2iT(Nf^ + 1). The size of the symbols codes for 
the magnitude of the quantity plotted and their position marks their x and y coordinates. Results for wave vectors with same x and y 
coordinates but differing z coordinates are overplotted on top of each other. Panels (a) and (b) compare the simulations with cubic cells 
(a) to the simulation with elongated cells (b) in the case where N]^=lb. Panels (c) and (d) compare numerical viscosity (c) and resistivity 
(d) in the case with elongated grid cells for N}^ = 5. 



8.4 Total pressure gradient 

The assumption of a homogeneous total pressure is the cornerstone of our analysis. This requirement might not be as strong 
as it seems at first glance. Indeed, the gradient of total pressure in the Euler equations naturally drives MHD flows towards 
a state of uniform total pressure. As a result, our solutions should be close approximations to the exact MHD flows even in 
cases when the condition of homogeneity is not met, provided that the real flow remains at a nearly constant density. 



9 CONCLUSIONS 

This paper consists of variations on a theme: the channel solution. We have extended previously known analytical solutions 
to more general and more physical cases, including viscosity, resistivity and cooling. We also showed the connection between 
torsional Alfven waves and channel solutions. 

We used these solutions to calibrate the implementation of a conservative scheme in Zeus3D. We also measured the 
numerical resistivity and viscosity of torsional Alfven waves in ZeusSD. In particular, we showed that lower time steps should 
be used in ZeusSD in order to guarantee a positive resistivity when no physical resistivity is used. We would rather recommend 
to use a minimal amount of physical resistivity. It also is best to use isotropic resolution since the numerical dissipation is more 
anisotropic for elongated cells. Finally we find a dependence of the numerical dissipation on the amplitude of the magnetic 
field. 

Although in this paper we stressed the numerical applications of these solutions, they are of interest in their own right. 
In particular, we have established a stronger basis for understanding the stability analysis of channel solutions: it should 
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now be possible to compute parasitic instabilities with improved microphysics. As a result, we hope to better understand the 
saturation properties of MRI turbulence. 
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APPENDIX A: PROPAGATING DISTURBANCES WITH ZERO SHEAR 



We work here in the framework and notations of section IT2l We first assume that there exists a common complex root s to 
R and Q. Then its complex conjugate s also is a common root. Being only second degree, R and Q have the same roots s and 
s. Hence they differ only by a real proportionality constant, and the remainder Ri of the Euclidian division of -R by Q needs 
to be identically zero: 

Ri{s)=ais + ao (Al) 
with 

ai = iy~v-7^ (A2) 



2 



and 



ao = {jy- v)v + k^vl 1^1 + ^^^"'' J . (A3) 

Both coefficients ao and ai must vanish if R and Q are to have a common complex root. This puts two constraints on the 
three remaining parameters i], v and kvA so that r\ and v can be expressed in terms of kvA- Setting a\ equal to zero, we get 

(A4) 

We now use this expression into the equation ao = which yields a quadratic equation for v: 

14!/^-2fc\ia!/ + 15A:\i =0 (A5) 
which has only one real positive root 



-a + 277 



1 

14 



— \ci^ Jcfi ^2\(Sk-^v\] . (A6) 



Equation HA4|) now provides the value for r\ 



r? = ^ ( 6a - ^q2 + 210fc2u^ ) (A7) 



which is positive for a > \f%\kvA\ ■ 

In order to get a common complex root to P and Q, the resistivity and viscosity need hence to be determined by 
expressions (IA7[) and (IA6[l . Using these expressions for r\ and v in the dispersion relation i?(s) = provides the two actual 
growth rates s and s as 

-17Q-3Va^ + 210fc2„2 

s± = ^ ± (A8) 

with 



r = 95k^VA + a^a^ + 210Pw^ - a. (A9) 

The corresponding solutions are therefore propagating disturbances (i.e. they have a non zero imaginary part) only when 
r > which is equivalent to setting the condition a > ~^\kvA\- Since we already required the more stringent condition 
a > ^/&\kvA\ in order to get 77 > 0, physically plausible solutions exist only when a > y/E\kvA\ for propagating disturbances. 
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APPENDIX B: STANDING WAVES WITH ZERO SHEAR 

We work here in the framework and notations of section 14.21 We now assume that s is a real common root to P and R. In 
that case, Ri{s) — immediately gives s in terms of the parameters of the other problem; 

_ -2iyri{u ~rj)+ a^vl - 2{r, + u)k\l 

Using s in P or Q, we finally arrive at a relational constraint for the defining of the problem. As an explicit example, k'^v\ 
can be expressed in terms of a, rj and v. 

k^v\ = l^-Q^ + 7Q(r7 + iy)+ Iv^ - 41?)!/ - 10?)^ + (-a + 2;/ + ^r])^J oP- - a(Vdv + 4?7) + + 44771/ + 4772 j . (B2) 



APPENDIX C: A METHOD TO MEASURE NUMERICAL RESISTIVITY AND VISCOSITY 



In section 14.11 we showed that circularly polarised waves with nonzero viscosity or resistivity are solution of the non-linear 
equations. In principle, if we start our simulation with one of the eigenmodes corresponding to the growth rate s± = ±iw 
(equation |4l] for ri = 1/ = 0), we should obtain the time evolution of a torsional Alfven wave as a result of the computation. 

However, the finite grid and time stepping resolution introduce some numerical defects. For example, figures [1^ and[T]D 
show that numerical results undergo some dissipation. In these figures, the dashed line corresponds to a wave with a slightly 
lower amplitude than a pure torsional Alfven wave after three oscillation periods. This suggests that the numerical errors in 
the code may behave like an equivalent viscosity and resistivity. In principle we could define their effective values if we were 
able to fit a model evolution to the actual numerical output of the code. 

In this appendix, we are motivated to compute the evolution of a system which starts with the initial conditions for a 
torsional Alfven wave (with rj — 1/ = 0), but which is evolved with some amount of viscosity z^n and resistivity t^n- Recall that 
77 — A;^77n and u = fc^i/N where 77N and are the effective resistivity and viscosity. We present the results to first order in 
(77 - v)/kvA = iVN ~ vn)k/vA- 

We choose the initial phase such that Sbx ~ 1. The initial conditions for a torsional Alfven wave give Sux = 1 and 
6by = Suy = i. 

We first assume that the code preserves well the initial uniform density profiltQ. According to section \4A\ there exist only 
two incompressible modes that can be excited with growth rates given by equation (|41|l . We decompose our initial conditions 
on the two corresponding eigenmodes which have 5ux± = (s± + ri)/{ikvA) Sbx±- 

oUx = 1 = a+ — ; h a- — : (^1) 

iw iw 

and 

56^ = l = Q++a- (C2) 
with 

^^^±1 (C3) 

and where a+ and q_ are the complex weights of the two eigen modes. 
We solve for a+ and q_ and retain the first order in (77 — u)/kvA'- 

1 / kvA+i^ \ , .ri-v 



and 



a„ = l-a+~-7|^. (C5) 

The non-linear coupling between these two modes can only occur through the total pressure gradient term and it happens 
that this term vanishes to first order in (77— u)/kvA- The temporal evolution of the system can hence be approximated by its 
linear evolution 

Sux = Q+ ^ exp(s+t) -I- a- - — ^tiZ exp(s_t) (C6) 
ikvA ikvA 

and 

^ We actually checked that to enforce p = 1 in the code did not change much the measured v and rj. 
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Figure CI. On the left panel (a) we display (solid line) the evolution of the contrast < (5)i[M])^ - (5R[i>])^ > / < (3?['"])^ + (3?[J>])^ > in 
our standard run as well as (dashed line) the result of equation l)Cll|l where rj — u is determined from the first time step of the simulation. 
On the right panel (b) we display (solid line) the evolution of the quantity 1— < (3?[ix])^ -I- (3?[b])^ > / < 2cos^(fc.r) > in our standard 
run as well as (dashed line) the result of equation l|C10|l where r] + u is determined from the first time step of the simulation. 



Sbx = a+ exp{s+t) + exp(s_f). 

We finally recover the temporal evolution of the perturbed quantities as 



and 



exp 



5R[M =exp(-^i 



cos{kvAt + k.r) 



cos{kvAt + k.r) 



2kvA 



sin{kvAt) cos(fc.r) 



2kvA 



sm{kvAt) cos(fc. 



The y component of these fields can be recovered because of the circular polarisation conditions by = ib^ and Uy = iu, 
choose to recover t) and u from their sum and difference through the volumic averages of kinetic and magnetic energy: 

\ < mu]f + >= cxp (-(r; -t- u)t) . 

and 

sin{2kvAt) 



(C7) 



(C8) 



(C9) 



We 



(CIO) 



2kvA 



< {^[u]y + {^[b]y > 



(CU) 



Hence, we find the total dissipation rj + u from the exponential decay of the kinetic plus magnetic energy. And we get the 
difference rj — v from the relative difference between these two forms of energy. Note that these final expressions yield rj + v 
with one more order of accuracy in {rj — v)/kvA than rj — u. We checked a posteriori that (77 — v)/kvA is indeed small for all 
measurements performed in this paper except for the highest wave numbers (see section [7T} . 

In order to save computing time, we evaluate rj and v on the very first time step of the simulation. Figures IClb , and IClb 
show that this provides reasonable estimates for equations (ICllI) and (|C10[) . The match is in fact perfect only for the first 
few time steps. The later discrepancy between our model for the code diffusion and the actual results of the code is probably 
due to the dispersive properties of the scheme which we do not take into account. This discrepancy actually narrows down at 
higher Courant numbers which are known to be less dispersive. However, figure CI shows that our model captures the bulk 
of the numerical artifacts. 

Finally, in order to ensure that the initial fields have zero divergence, we convert real wave numbers k to discrete wave 
numbers k' — 2sm{kAx /2) / Ax when we compute the relations between the amplitudes of the fields. Ax is the size of a pixel 
in the units of the computation. 
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